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I .      INTRODUCTION 

A  basic   building  block  for  describing  underwater   explosions   is   the 
hydrodynamic   properties  across  a   shock  wave   in  sea  water.      Underwater 
explosions  were  very  extensively  studied   in  World  War   II.      Results   from 
the  World  War   II  studies  are  reported   in  Underwater   Zxplosion  Research    [1] . 
Richardson,    Arons,   and  Halverson    [2]   discuss    the  calculation  of  hydro- 
dynamic   properties  of   sea  water  at   the  front  of  a   shock  wave.      Background 
information  on  the  calculations  appears   in  Chapter   2  of  Cole   [3]. 

Richardson,    et  al. ,    [2]    used  graphical   techniques   to   interpolate 
thermodynamic  data;    the  methods  were  crude  and   extremely   tedious.      With   the 
advent  of   the  modern  programmable  hand-held   calculator,    the  results 
obtained  by  Richardson,    et  al.,    [2]    can  be  obtained  with  ease.      This  report 
discusses  a   program  for   the  HP41CV,   which  uses    the  data  of   Richardson, 
et  al. ,    [2]    to  calculate  shock  wave  velocity,   water  velocity,    specific 
volume,    etc.      Reference    [2]    is  Appendix  A. 

Since  World  War   II,    research  has  continued   in  the  area  of   underwater 
explosions.     A  more  recent  survey  by  Professor  Holt    [4]   appears   in  the 
Annual  Reviews  of   Fluid  Mechanics.      Holt    [4]   discusses  underwater   nuclear 
explosions  as  well  as  chemical  explosions. 

In  passing,   an  interesting  observation  can  be  made  concerning   the  level 
of   talent  working  on  underwater   explosions   in  World  War   II.      The  names 
include  G.    I.    Taylor,    P.    Bridgman,    H.    Bethe,    J.   von  Neumann,   and  L.    I.    Sedov 


II.   OUTLINE  OF  CALCULATION  PROCEDURE 

A.   Iteration  of  Properties  with  Rankine  Hugoniot  Relations 

The  computer  program  has  an  iteration  which  is  now  discussed.   The 
enthalpy  jump  across  the  shock  wave  is  given  by 

AH  =  |(p  -  pQ)(v  +  vQ)  (2.3) 

The  same  equation  numbers  and   the  same  symbols   used  by  Richardson,    et  al. ,    [2] 
are  used  here.      Since  p  >>  p_,   p     can  be  neglected.     Also  AH  is  given  by 

AH  -  a)  +  h  (2.20) 

where  oj  is   the  undissipated   enthalpy  and  h  is   the  dissipated   enthalpy.      The 
dissipated   enthalpy  can  be  calculated   two  ways.      First, 

hy  =»  AH  -  oj 
Second, 

Tl 
AH  -  h     -  //  C  dT  (2.9) 

P  T         TQ     p 

The  symbol  tv,,   which  does  not  appear   in  reference   [2],   denotes   the  value 
calculated  from  AH  -  u).     The  symbol  h     is   the  value  calculated   from 
equation   (2.9).      Temperatures  TQ  and   T.,   are  defined   in  reference    [2].      The 
undissipated  enthalpy  is  a  function  of  specific  volume,   v,   behind   the  shock 

wave;    the  function  is 

2  n-1 

c  v 

u  -T^rt  [CT°     ~  1]  (2-19) 

The  Tait  equation  of   state  relating  specific  volume  and  pressure  for  water   is 

P  -  AtSlUv^v)11-  1] 
where  A[S]    is  a  function  only  of   entropy.      Further 

A[S]    =•  c^/nv1  (2.14) 


Combining  equations  (2.3),  (2.19),  and  (2.14)  yields 

2  2 

c  ~  c  ~ 

h  =  2n7~  [(vi/v)tl  "  1](v  +  V    "  n~^tC(vl/v)n"1   "  1]  (2>21) 

Equation   (2.21)    indicates   that  h_   is  a   function  of  v  and  v    .      The  specific 

-i  1 

volume,  v  ,  is  a  function  of  T  .   By  varying  T  and  v,  the  value  of  h^ 
varies. 

When  conditions  across  the  shock  wave  are  correct, 

hT(T1)  =  hH(v,T1) 

One  varies  v  and  T  until  equality  is  achieved;  since  v  =  v(p,T  ),  only  T- 
needs  to  be  varied. 

This  is  a  very  brief  description  of  the  iteration.  The  reader  should 
refer  to  reference  [2]  for  more  information. 
B.   Temperature  Behind  Shock  Wave 

The  temperature  increase  behind  the  shock  wave  was  calculated  using 
equation  (II-4)  of  reference  [2] .   Although  the  equation  is  somewhat 
lengthy,  the  calculation  is  straightforward  except  for  two  quantities: 
B'Ctfl)  and  8n-  A  cubic  polynomial  is  given  for  3(tn);  this  polynomial  was 
differentiated  to  obtain 

B'(tJ  -&) 

u    dt 

t=t0 

The  quantity   3n   is    the  coefficient  of    thermal   expansion  times   specific 


volume  and  was   evaluated   using 


Av 
S0   "It 


where  v  is  specific  volume.   The  symbol  t  is  used  for  C  and  T  for  K. 


III.   DESCRIPTION  OF  COMPUTER  PROGRAM 

A.   Computer  Program  Plow  Chart 

A  flow  chart  of  the  program  is  given  in  Figure  1.   An  iteration  loop 
occurs  for  making  h„  3  HH  =  h_  ■  HT.   The  input  "LOOP"  is  the  criterion 
for  acceptable  difference  between  HH  and  HT.   Since  HH  and  HT  have  units 
of  Joule/kg,  "LOOP"  also  has  units.   The  iteration  variable  is  t  ,  which 
is  the  adiabatic  temperature  in  C. 

When  an  acceptable  value  for  t.  has  been  found,  the  program  auto- 

3 
matically  goes  to  subroutine  SHOCK.   SHOCK  calculates  v,  cm  /gm;  u,  m/sec; 

U,  m/sec;  and  c,  m/sec. 

Following  SHOCK,  the  program  moves  to  subroutine  DEL  T,  which  calculates 
the  difference  in  temperature  from  behind  to  front  of  shock  wave.   A  variety 
of  subroutines  are  used  by  DEL  T  including  CP,  BTPRIME,  Gl,  and  D2. 
3.  Assignment  of  Storage  Registers 

For  the  program,  34  storage  registers  are  used.   The  assignment 
of  variables  to  the  registers  is  given  in  Table  I. 
C.   Program  Listing 

Appendix  B  is  a  listing  of  the  program. 

IV.   OPERATION  OF  PROGRAM 

The  program  has  been  recorded  on  magnetic  cards  and  can  be  inserted 
into  the  HP41CV  using  a  HP82104A  card  reader.  The  inputs  to  the  program 
will  now  be  described. 

1.  Turn  on  the  HP41CV  and  activate  USER. 

2.  Assign  the  program  to  a  key  on  the  key  board  by  pressing 

□  ALPHA  HH  ALPHA,  SIN 


INPUT 

LOOP 

INPUT 

N=n 

YES 


GO  TO  "SHOCK" 


CALCULATE 
a+1  .   1    1 
n-1      n    n 


CALCULATE 

SPECIFIC 

VOLUME  v  BEHIND 

SHOCK 

WAVE 

INPUT  INITIAL 
TEMPERATURE 


IN°C 


INPUT  FINAL 
PRESSURE  IN  KILOBARS 


CALCULATE  v 
XEQ  VW 


'! 


INPUT  t. 


XEQ  Y 
CALCULATE  (v  /v)° 

FOR  USE  IN  EQ.  (2.21), 


NO 


XEQ  VW 
CALCULATE  v, 


XEQ  SPEED 
CALCULATE  c^/n 


XEQ  HT 
CALCULATE  h. 


CALCULATE  h_ 

DISPLAY  HH,  HT 
AND  HH-HT  = 

1 

ABS (HH-HT)  <  LOOP 

ir 

YES 


CALCULATE  SPEED 

OF  WATER  BEHIND 

SHOCK  WAVE,  u 


CALCULATE  SPEED 
OF  SHOCK  WAVE,  U 


CALCULATE  SPEED 

OF  SOUND  3EHIND 

SHOCK  WAVE,  c 

!  XEQ  DEL  T  ! 


XEQ  CP 
CALCULATE  C 


XEQ  BT  PRIME 
CALCULATE  dB/dt 


XEQ  Gl 

CALCULATE  G 


XEQ  D2 
CALCULATE  D 


CALCULATE  AT 


DISPLAY  AT 


END 


Figure  1.   Computer  Flow  Chart. 


"able  I.   Storage  Registers 


Register   Symbol 


Definition 


Units 

Programs 

°C 

HT,  VW,  HH, 
SPEED 

°C 

3T 

°c 

VW 

°c 

BT 

°c 

HH 

Joul< 

s/kg 

HT 

Joul< 

a/kg 

HH 

01 


19 


adiabatic  temperature 


02 

(t-55r 

quantity 

in  polynomial 

03 

(t-25) 

quantity 

in  polynomial 

04 

(t-55) 

quantity 

in  polynomial 

05 

C0 

temperature  in  front  of 
shock  wave 

06 

hT 

HT,  dissipated  enthalpy 

07 

h 

HE,  dissipated  enthalpy 

08 

B(t) 

quantity 
of  state 

in  Tait  equation 

09 

3(t) 

quantity 
of  state 

in  Tait  equation 

10 

t-55 

quantity 

in  polynomial 

11 

A[S] 

quantity 
of  state 

in  Tait  equation 

12 

P 

pressure 

behind  shock  wave 

13 

V 

specific 

volume 

14 

u 

velocity 

of  water  behind 

shock  wave 

15 

p 

pressure 

behind  shock  wave 

16 

U 

velocity 

of  shock  wave 

17 

Cl 

speed  of 

sound  at  adiabatic 

temperature 

18 

c 

speed  of 
wave 

sound  behind  shock 

coefficient  of  thermal 
expansion  times  specific 
volume  v„ 


kilobars   BT 


N/m 
°C 

N/m' 


N/m 
m/sec 

m/sec 

m/sec 


BT 

BT  PRIME 

SPEED 


kilobars  Y 

cm3/gm  SHOCK,  VW 

m/sec  SHOCK 
2 


SHOCK 
SHOCK 

SHOCK 

SHOCK 


m  /kg°K    BETA 


Table  I  Continued.   Storage  Registers 


Register   Symbol 


Definition 


Units 


Programs 


20 

c./n 

quantitv   in  equation 
(2.21) 

m"/sec" 

SPEED 

21 

(n+l)/(n- 

-l) 

function  of  n 

- 

HH 

22 

1-1/n 

function  of   n 

- 

HH 

23 

-1/n 

function  of  n 

- 

EH 

24 

dB/dt 

derivative  of  B(t) 

N/m2  °K 

BT  PRIME 

25 

y 

(v./v)      ratio  of 

specific  volumes 

- 

HH,    Y 

26 

vi 

specific  volume  at 
adiabatic   temperature 

cm  /gm 

HH,    VW 

?7   v 
vQ 

28  n 

29  C 

P 

30  G 

31  D 

32  1+p/B 

33  (1+p/B) 

34  LOOP 


1/n 


specific  volume  in 
front  of  shock  wave 

exponent  in  Tait  equa- 
tion of  state 

heat  capacity  at 
constant  pressure 

quantity  in  equation 
(II-4) 

quantity  in  equation 
(II-4) 

quantity  in  Tait  equa- 
tion of  state 

quantity  in  Tait  equa- 
tion of  state 

criterion  for  acceptable 
difference  between  HH 
and  HT;  typical  value 
1.0 


cm  /gm 


Joule/kg 


HH 

HH 

C? 

G2,  DEL  T 

D2,  DEL  T 

HH 

HH 


Joule/kg   HH 


In   this   example,    the  program  has  been  assigned    to    the  SIN  key 
with  blue  H. 

3.  Press   SIN  and  observe  LOOP?.      LOOP     is    the  value  of  HH-HT  which 
is  acceptable.     When  ABS(HH-HT)    is   less    than  LOOP,    the  program 
is  out  of   the  iteration  loop.     A  typical  value  for  LOOP 

could  be  1.0  Joule/kg.      Values  for  a   sample  case  are  given  in 
parentheses.      (1.0) 

4.  Press  R/S  and  observe  N? .      N  is   the  exponent  in  Tait   equation 

of  state.      Holt    [4]   uses   7.0.      Richardson,    et_al . ,    [2]    use  7.15. 
Insert  a  value  for  N.      (7.15) 

5.  Press  R/S  and  observe  TEMP  0?.      This   is   the  temperature  in  front 
of   the  shock  wave  in     C.      Insert  a  value.      (0) 

6.  Press  R/S  and  observe  PRESSURE?.      This   is   the  pressure  behind  the 
shock  wave   in  kilobars.      Insert  a  value.      (80) 

7.  Press  R/S  and  observe  TEMP   1?.    This   is   the  adiabatic   temperature 
discussed   in  reference   [2],      Insert  a  value  in     C.      (180) 

8.  Press  R/S  and  observe  crows  foot  moving  across   the  display.      In 
10  seconds,   a  value  for  HH  appears.      (520,245) 

9.  Press  R/S  and  observe  value  for  HT.      (724,255) 

10.  Press  R/S  and  observe  value  for  HH-HT.      (-204,010) 
Since  the  value  for  HH-HT  exceeds  LOOP,   a  new  value  for 
t     3  TEMP  1  is  needed. 

11.  Press  R/S  and  observe  TEMP  1?.      To  decide  what  value  of   TEMP  1 
to  insert,   look  at  equation  (2.9);   h     increases  as  T.    or  TEMP   1 


increases.      Although  not   readily  apparent   from  equation   (2.21), 
h_  decreases   as   TEMP   1   increases.      For    the   example,    a  small 
table  can  be  made  as   follows: 


TEMP  1 

HH 

HT 

HH-HT 

180 

520,245 

724,255 

-  204,010 

160 

569,362 

642,346 

-  72,984 

149.2 

591,574 

598,988 

-  7,414 

148.0 

593,821 

594,120 

-  299 

147.75 

594,285 

593,106 

1,180 

147.949 

593,916 

593,913 

3 

147.9496 

593,915 

593,915 

-  0 

12.  Press   R/S  and  observe   the  following: 

3 
v  *  0.6576      (cm  /gm) 

WATER  U  =   1,638      (m/sec) 

SHOCK  U  =  4,891      (m/sec) 

WAVE  C  -  6,264      (m/sec) 

13.  Press  R/S  and  observe   the  following: 

DEL  T  »   371.24°K 
Hence  the  temperature  behind   the  shock  wave  is 
T  ■  TQ  +  AT  -  273.16  +  371.24  =  644. 4°K 
The  calculation  of  shock  properties  has   now  been  completed.      To  calculate 
values   for  another   set  of    t_  and  p,    press   SIN  key  with  blue  H.      The 
computer  is   in  USER  mode. 


V.  SAMPLE  RESULTS 

A  comparison  can  be  made  between  results  calculated  with  the  program 
and  values  tabulated  by  Richardson,  et  al. ,  [2].   Table  II  is  a  comparison. 
Most  values  are  within  1/2  per  cent. 

Since  Richardson,  et  al. ,  [2]  used  graphical  techniques  for  interpola- 
tion, the  calculator  values  probably  are  more  accurate.   However,  the 
thermodynamic  data  for  water  have  been  represented  by  polynomials.   Errors 
in  curve  fitting  are  undoubtedly  a  few  per  cent.  Hence,  one  would  expect 
the  values  predicted  here  to  be  correct  within  a  few  per  cent. 

VI.  MAGNETIC  CARDS 

Copies  of   the  magnetic  cards  containing   this  program  are  available 

on  request  from  the   following  address: 

Distinguished  Professor  Allen  E.    Fuhs 
Department  of  Aeronautics,    Code  67 
Naval  Postgraduate   School 
Monterey,    CA  93940 

(408) -646-2943 
AV-8 78-2948 

VII.   CALCULATION  OF  WATER  PROPERTIES 
DURING  ISENTROPIC  EXPANSION  FROM  SHOCKED  CONDITIONS 

In  order  to  calculate  the  complete  flow  field  due  to  a  shaped  charge 
jet  penetrating  water,  the  calculation  of  temperature,  internal  energy,  and 
enthalpy  during  isentropic  expansion  is  necessary.   Appendix  C  discusses  two 
subroutines  which  were  developed  after  the  main  program  was  written.  The 
subroutines  calculate  internal  energy  and  enthalpy,  subroutine  E-3AR,  and 
temperature,  subroutine  ISEN  T. 
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APPENDIX  k 
HYDRODYNAMIC  PROPERTIES  OF  SEA  WATER  AT  THE  FRONT  OF  A  SHOCK  WAVE* 


*This  is  reference  [2]  reproduced  from  reference  [1].   The  technical  paper 
is  provided  as  a  convenience  to  the  reader. 
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Hydrodynamic  Properties  of  Sea  Water  at  the  Front  of  a  Shock  Wave* 

J.  M.  Richardson** 

Baker  Laboratory,  Cornell  University,  Ithaca,  .V.   Y. 

AND 

A.  B.  Aeons***  .and  R.  R.  Halverson**** 
Underwater  Explosives  Research  Laboratory,  Woods  Hole  Oceanographic  Institution,  Woods  Hole,  Massachusetts 

(Received  June  2,  1947) 

The  Rankine-Hugontot  relations  have  been  applied  to  appropriate  equation-of-state  data 
in  order  to  calculate  the  propagation  velocity,  particle  velocity,  enthalpy  increment,  Riemann 
function,  etc.  at  shock  fronts  of  various  amplitudes  in  sea  water.  One  set  of  tables  provides 
values  over  a  wide  pressure  range  vup  to  about  80  kilobars)  and  is  principally  intended  for  use 
in  conjunction  with  theories  of  propagation  oi  shock  waves  originated  by  underwater  ex- 
plosions. A  second  set  of  tables  contains  values  which  are  closely  spaced  up  to  pressures  of  14 
kilobars.  These  are  calculated  with  somewhat  greater  precision  and  are  intended  for  use  in  con- 
nection with  experimental  measurements  of  particle  and  propagation  velocities,  etc. 


L  HTTRODUCTIOIT 

IT  has  long  been  recognized  that  the  velocity 
of  propagation  of  sound  waves  of  finite 
amplitude  in  a  fluid  medium  is  a  function  of  the 
pressure  in  the  wave.  Lamb1  ascribes  the  early 

•  The  work  described  in  this  report  was  performed 
under  National  Defense  Research  Committee  Contracts 
OEMsr-121  with  Cornell  University  and  OEMsr-569  with 
the  Woods  Hole  Oceanographic  Institution. 

*•  Present  address:  Bell  Telephone  Laboratories.  Inc.. 
Murray  Hill.  N.  J. 

""Present  address:  Department  of  Physics.  Stevens 
Institute  of  Technology,  Hoboken,  N.  J. 

•**•  Deceased. 

1  H.  Lamb,  Hydrodynamics  (Cambridge  University  Press. 
London,  1932)  6th  Ed.,  p.  481. 


development  of  the  theory  to  independent  in 
tigations  of  Earnshaw  and  Riemann.  Qua. 
tively  this  work  indicated  that,  since  the  hij 
pressure  portions  of  a  wave  travel  with  gre 
velocity,  an  arbitrarily-shaped  pressure  puis 
finite  amplitude  must,  during  propagation,  a 
its  shape  in  such  a  manner  as  to  build  up  in 
shock  front.  By  applying  the  laws  of  conse: 
tion  of  mass,  energy,  and  momentum  to 
transfer  of  matter  across  the  shock  fr 
Rankine  and  Hugoniot  obtained  a  set  of  tl 
relations  among  the  five  variables:  press 
density,     particle     velocity     («),     shock     fi 
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Fig.  1.  Adiabatic  and  Hugoniot  contours  in  p—  T  plane. 

velocity  (£/),  and  enthalpy  increment  (AH). 
These  relations,  when  applied  to  data  on  the 
equation-of-state  and  specific  heat,  make  it  pos- 
sible to  calculate  u,  U,  and  iH  and  to  evaluate 
certain  other  functions  applicable  to  the  theory 
of  the  formation  and  propagation  of  shock  waves 
originated  by  explosions.1 

Precise  knowledge  of  u  and  U  also  makes  it 
possible  to  calculate  shock  wave  pressures  in 
cases  where  the  particle  velocity'  or  propagation 
velocity  can  be  measured.  The  purpose  of  the 
calculations  described  betow  was  to  apply  the 
Hugoniot  relations  to  appropriate  equation-of- 
state  data  for  sea  water  in  order  to  provide  a) 
tables  of  the  desired  functions  up  to  very  high 
pressures  {ca.  80  kilobars)  for  use  in  the  theory 
of  propagation  of  underwater  explosion  waves,2 
and  (b)  tables  of  particle  and  propagation 
velocity  at  fairly  close  pressure  intervals  in  a 
lower  pressure  region  (up  to  ca.  14  kilobars). 

DL  OUTLINE  OF  THE  THEORY  AND  COM- 
PUTATIONAL PROCEDURES 

In  this  section  we  give  an  account  of  the 
hydrodynamical  and  thermodynamical  relations, 
and  the  computational  procedures  leading  to  the 
numerical  results  tabulated  in  Sections  III  and 
IV.  For  the  convenience  of  the  reader  a  glossary' 
of  symbols  is  presented  in  Appendix  III. 

When  a  shock  wave  advances  with  velocity  U 
into  a  stationary  fluid  of  unperturbed  pressure 

*  J.  G.  Kirkwood  and  H.  Bethe.  The  Pressure  Wave 
Produced  by  an  Underwater  Explosion  Dept.  of  Commerce 
Bibliography  No.  PB  J2182),  OSRD  Reoort  No.  588. 
Part  I. 


po  and  specific  volume  :-0,  the  pressure  p.  specific 
volume  ;•,  and  particle  velocity  u  of  the  rluid 
behind  the  shock  front  are  determined  by  the 
Rankine*- Hugoniot4  conditions,  which  express 
the  conservation  of  mass,  momentum,  and 
energy  of  an  element  of  fiuid  passing  through  the 
front.  For  the  purposes  of  this  paper,  these  con- 
ditions mav  convenientlv  be  written 


«  =  [(£— pt)(vo—  «)]*, 

U  =  volip  -  p0)  /  (»o  -*)]*■ 

\H=n/2)(p-po)(v^v<,). 


(2.1) 

(2.2) 
(2.3) 


In  the  last  equation,  AH  is  the  specific  enthalpy- 
increment  of  an  element  of  rluid  when  it  passes 
through  the  front.  The  specific  enthalpy  is  de- 
fined as  the  sum  of  the  internal  energy  per  gram 
and  the  pressure-volume  product,  pv. 

Given  equation  of  state  and  specific  heat  data 
for  the  fluid,  any  three  of  the  variables  p.  v,  U 
and  u  may  be  determined  as  functions  of  the 
fourth.  Here  we  shall  regard  p  as  the  independent 
variable.  For  certain  hydrodynamic  applications 
we  must  have,  in  addition  to  v,  U,  and  u  as 
functions  of  p,  the  sound  velocity 


c  =  ■  dp/dp)*  s;     p=l,v, 
the  Riemann  <r-functionf 

ir»  f    [?~J>' .  S'},  c\_p' ,  S7]dp' , 

and  the  undissipated  enthalpy 

u-  f  v\j',syLp\ 


i,2.4) 


!.5) 


(2.6) 


where  5  is  the  entropy. 

In  practice,  one  must  resort  to  successive 
approximations  to  effect  a  reduction  of  the 
Hugoniot  conditions,  combined  with  equation-of- 
state  and  specific  heat  data,  to  a  set  of  relations 
expressing  u,  U,  and  v  as  functions  of  p.  To  this 


Roy.  Soc.  London,  A160. 

1387);  58,   ! 


1  W.  J.   M.  Rankine.  Trans 
277  (1870). 

•  H.  Hugoniot.  J.  de  l'ecole  polvt.  51 
1888!. 

+  The  Riemann  ir-tunction  occurs  in  Riemann's  form  of 
the  hydrodvnamicai  equations,  which,  for  the  case  oi 
spherical  symmetry,  may  be  written    see  reference  I): 


i'd/dt)-r-<c  +  u)(a/dr)V<r-r-u)-- 
Tjdidt)  —  (c  —  uHd/dr)  *(<r—  u)> 


<  —  Zcu/r, 
•0. 


where  i  is  the  time  and  r  is  the  radial  coordinate;  The 
other  quantities  have  already  been  defined. 
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end,  it  is  expedient  nrst  to  consider  certain 
quantities  as  functions  of  pressure  and  tem- 
perature, p  and  T,  or  pressure  and  entropy,  p  and 
J.  Before  proceeding  to  a  more  detailed  dis- 
cussion of  the  calculations  it  may  perhaps  heip 
to  orient  the  reader  if  we  consider,  qualitatively, 
contours  of  some  pertinent  quantities  in  the  p-T 
plane.  In  Fig.  1,  the  possible  states  of  a  given 
fluid  just  behind  the  shock  front  lie  along  a  single 
curve,  which  we  have  labeled  "Hugoniot."  An 
element  of  fluid  initially  in  the  state  lpo,  To) 
which  has  attained  a  state  (J>,  T)  just  behind 
the  shock  front  finally  returns  to  a  state  {pa,  I\) 
along  the  adiabatic,  so  labeled  in  the  figure.  Also 
included  are  the  designations  of  a  few  points  on 
a  p-S  basis  using  square  brackets  according  to 
the  convention  introduced  in  Part  b  of  this 
section.  In  general,  7\  is  larger  than  To  because 
of  the  dissipation  occuring  at  the  front.  The 
central  part  of  our  problem  is  the  determination 
of  the  Hugoniot  curve. 

We  shall  consider  in  Part  a  the  calculations 
due  to  Arons  and  Haiverson*  which  are  intended 
to  be  accurate  in  the  range  of  relatively  low 
pressure  (ca.  0  to  20  kilobars).  These  results  as 
stated  in  the  introduction  are  intended  for  the 
determination  of  the  peak  pressure  of  a  shock 
wave  from  measured  values  of  the  shock  front 
velocity  U  or  particle  velocity  u.  In  Part  b,  we 
shall  consider  the  calculations  of  Kirkwood  and 
Richardson,4  the  results  of  which  were  originally 
intended  for  the  applications  of  the  shock  wave 
propagation  theory  of  Kirkwood  and  Bethe* 
which  required  data  over  a  higher  range  of  pres- 
sures (ca.  20  to  50  kilobars). 

a.  Calculations  of  Arons  and  Haiverson 

Here  we  outline  the  calculations7  suitable  for 
the  relatively  low  pressure  range  [ca.  0  to  20 
kilobars)  based  upon  the  equation-of-state  and 
specific  heat  data  discussed  in  detail  in  Appendix 
I.  For  the  range  0  to  1.5  kilobars,  the  Ekman 
equation-of-state  was  used ;  in  the  range  0  to  25 

*  A.  B.  Arons  and  R.  R.  Haiverson,  Hueoniot  Calcula- 
tions for  Sea  Water  at  '.he  Shock  Front,  OSRD  Report 
No.  6577,  NDRC  No.  A--M59. 

'J.  G.  Kirkwood  and  J.  M.  Richardson,  The  Pressure 
Wave  Produced  by  an  Underwater  Explosion,  Part  III,  OSRD 
Report  No.  813  I  Dept.  of  Commerce  Bibliography  No.  ?B 
32184). 

TJ.  G.  Kiricwood  and  E.  M  on  troll.  Pressure  Wave 
Produced  by  an  Underwater  Explosion,  II,  OSRD  Report 
No.  676  vDept.  of  Commerce  Bibliography  PB-32183). 


kilobars,  the  Tait  equation-of-state, 

v[0,  T)-v(p,  Z*)/f(0,  T)  - (l/«)  logCl-H>/5 

*  =  iT-273.16)°C. 

In  the  first  case,  the  initial  temperature 
fo=15°C;  in  the  second,  .'o  =  25"C.  In  both  c 
the  initial  pressure  po=*0.  Neither  of  the 
equations-of-state  are  complete  in  the  sense 
i'o  =  f(0,  T)  must  be  determined  by  auxi 
thermal  expansion  data  'also  discussed  in 
pendix  I). 

We  express  the  enthalpy  and  volume  ii 
ments 

AH=H(p,  T)-H{Q,  To), 
±v  =  v(p,  T)-v{0,  To), 

in  terms  of  line  integrals,  first  along  an  is 
from  (.0,  To)  to  (0,  T)  and,  secondly,  aloni 
isotherm  from  (0,  T)  to  {p,  T)  (see  Fig.  1). 
the  enthalpy  increment  we  obtain 


J  r, 


Cp(0,  DdT'djAT, 


dv(p',  T) 

v(p',  T)-T 

l  dT 


dp', 


±T=T-Tu 


where  c„(0,  T)  is  the  specinc  heat  extrapolat« 
zero  pressure  and  <?„  is  the  mean  of  cv  ovei 
temperature  range  \T. 

For  the  volume  increment  we  obtain 

■V=r(0,  D  -v(0,  To)  =*$*iT,       (j 
±tv=*v(P,  I)—  0(0,  D, 

where  Jo  is  the  mean  thermal  expansion  at 
pressure  over  the  temperature  range  AT. 

From  the  last  Hugoniot  condition,  Eq.  ( 
and  Eq.  (2.9)  we  obtain 

C»(0,  r0)  +  (l/2)(Art;)>-Ar^ 

Ar- ,  a 

5,-(l/2)0V) 

where  ItH  is  to  be  calculated  by  means  of 
third  of  Eq.  (2.9)  and  the  appropriate  equal 
of-state,  and  where  At»  is  to  be  obtained  f 
compressibility  data.  The  right-hand  side  of 
(2.11)  depends,  of  course,  on  the  temperatur 
The  determination  of  Ar  is  accomplished  by 
method   of   successive  approximations.   A 
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Table  !.    L'—r9 

based   on    Ekman 

temperature   lS'C; 

wt.  percent  XaC!) ; 


.).  u.  md  \H  in  low  pressure  -e^ion 
;quauon-ot'-*rate'.  Sea  .vater:  initial 
-alinitv  32  cans  per  thousand  (3.79 
-.,  =  4922.3  ft,  ^ec=»  1500  5  m  sec. 


.4 

a 

V  —ca 

u 

A/7 

I 

V  —  c» 

* 

iif 

r-Pt 

£a 

•  m/ 

cal, 

•>  —  i>i 

:a 

(ft, 

caj  / 

Icbar* 

!%) 

secj 

gmi 

lb  in') 

(%) 

sec) 

?m. 

o.oo 

000 

0.0 

00 

0 

(.00 

0  0 

0.0 

.25 

:.o" 

15.1 

5.3 

2.J00 

l.U 

29  1 

3  2 

so 

4.03 

31.2 

11.5 

+.000 

2.27 

57.6 

6.4 

.75 

5.93 

46.2 

17.2 

6.000 

3JS 

J5.5 

96 

1  on 

7.81 

60.5 

23.0 

3.000 

+  .43 

112.9 

12.7 

1.25 

■J. 05 

7*.+ 

23.5 

10.000 

5.48 

139.7 

15.9 

1.50 

11.44 

S7.7 

34.1 

12.000 

0J3 

166.0 

19.0 

U.000 

7.56 

191.8 

22.1 

•■     16.000 

i.58 

217.1 

25.2 

18.000 

9  59 

2+2.0 

23J 

20.000 

10.57 

266.5 

31.4 

1    22.000 

11.55 

290.6 

3+5 

Table  II.  U—  coiCa.  u,  and  SH  in  intermediate  pressure 
.egion  1.5  to  25  kilobars)  based  on  Tait  equation-ol- 
-tate;  t  =  7.300,  3 «•  3.012).  Sea  water:  Initial  temperature 
25°C;  salinitv  32  parts  per  thousand  '  3.79  wt.  percent 
N'aCl):  co-5014.7  ft/sec- 1528.5  m/sec. 


4 

B 

IS  — c« 

u 

.iff 

U-ti 

H 

±H 

f  -71 

-» 

(m/ 

cal/ 

P—  P* 

-• 

':t/ 

cal/ 

Icbar  ( 

C%) 

«C) 

gm) 

Ab/  in') 

'  '  o  i 

sec) 

jmi 

oo 

0.00 

0.0 

0.0 

0 

0  0 

0 

0.0 

1.0 

8.18 

59.2 

23  0 

20.000 

11.0 

258 

31.3 

1.5 

11.38 

3S.5 

34.2 

30.000 

15.3 

375 

46.0 

2.0 

15.39 

110.1 

45.3 

40.000 

20.5 

483 

02.0 

2.5 

18.71 

135.0 

56.3 

50.000 

24.7 

534 

77.0 

3.0 

21.38 

157.8 

67.2 

60.000 

28.7 

676 

92.0 

4.0 

27.34 

200.5 

^8.8 

TOOOO 

32.5 

■67 

107.0 

5.0 

33.39 

240.2 

110.1 

30.000 

36.0 

JS3 

121.0 

6.0 

38.59 

277.5 

131.2 

TOOOO 

39.4 

937 

13S.0 

8.0 

48.13 

346.1 

172.9 

100  000 

+2.3 

1014 

150.0 

10.0 

56.73 

408.9 

214.0 

120  000 

49.2 

1168 

178.0 

12.0 

64.70 

466  9 

2S4.3 

140.000 

55.2 

1312 

207.0 

14.0 

72.25 

510.9 

295.2      , 

InO.OOO 

00.8 

1445 

235.0 

25.0 

108. 40 

769.1 

51+.+ 

.30.000 

66.3 

1564 

263.0 

200.000 

71.5 

1662 

291.0 

value  of  AT  is  used  in  evaluating  the  right-hand 
side  giving  a  more  accurate  value  of  ±T  on  the 
left-hand  side,  and  the  process  is  repeated  until 
the  results  of  two  successive  steps  differ  by  a 
sufficiently  small  amount.  One  or  two  steps 
generally  suffice. 

We  have  thus  obtained  T  as  a  function  of  p 
along  the  Hugoniot  curve  :see  Fig.  1).  It  is  now 
possible  to  calculate  immediately  the  particie 
velocity  z*,  the  propagation  velocity  V,  and  the 
specific  volume  v  as  functions  of  p  behind  the 
shock  front.  The  results  using  the  Ekman 
equation-of-state  and  the  Tait  equation-of-state 
are  tabulated  in  Tables  I  and  II,  respectively,  of 
Section  III. 

b.  The  Calculations  of  Richardson  and  Kirkwood 

Here  we  outline  the  calculations4  intended  for 
the  applications  of  the  shock  wave  propagation 


theory  of  Kirkwood  and  Berne. :  These  are  based 
upon  the  equation-of-state  and  specific  heat 
data  discussed  in  detail  in  Appendix  II.  We 
use  a  modified  Tait  equation-of-state  connecting 
vip,  T)  and  r(0,  T)  to  be  discussed  below.  In 
most  respects,  the  data  is  made  to  tit  the  proper- 
ties of  an  aqueous  0.7  molai  XaCl  solution 
assumed  to  be  roughly  equivalent  to  sea  water 
of  salinity  s  =  32  parts  per  thousand  <see  Section 
1  of  Appendix  I). 

In  these  calculations  the  initial  pressure  po  is 
taken  to  be  zero,  and  several  different  initial 
temperatures  T0  are  used :  0°C,  20°C.  and  40°C. 

Before  indicating  the  precise  nature  of  the 
modification  oi  the  Tait  equation,  it  is  desirable 
to  mention  that  in  this  part  two  different  pairs 
of  independent  variables  will  be  used :  pressure 
and  temperature  (p,  T),  and  pressure  and 
entropy  7_p,  Sj.  Consequently,  in  order  to 
indicate  which  pair  are  used  in  a  function,  we  will 
use  parenthesis  to  indicate  the  first  pair  and 
square  brackets  to  indicate  the  second,  i.e. 
vip,  T)  and  v[p,  5]. 

The  modified  form  of  Tait  equation  introduced 
by  Kirkwood2- 5  is 

log(tri/»)»(l/»)  log "v  1  —  p/AfS"]),     (2.12) 

where 

v-vfp,  5]=r(p,  7Tj>,  5]),     yl  =  f[0,  S], 

see  Fig.  1)  n  is  an  empirical  constant,  and  the 
function  A [5]  is  related  to  the  function  Bit)  in 
the  original  isothermal  form  of  the  Tait  equa- 
tion, Eq.  (2.7),  as  follows, 

4C5]-5(C0,5]),     t-(r-273.16)°C.     (2.13) 

The  reasons  for  introducing  this  modification  of 
the  Tait  equation  are  at  least  twotold :  (1)  the 
anomaly  of  a  vanishing  specific  volume  v{p,  T)  at 
a  finite  pressure  along  a  jiven  adiabatic  (which 
does  not  differ  markedly  from  the  Hugoniot  curve 
in  the  case  of  water)  is  removed  to  a  higher  pres- 
sure  by   replacing   [v{0,  T)  —  vip,  T)~]/[v{$,  T)~\ 

Table  III.  Values  ot  U—co/c*  tor  different  temperatures 
and  salinities  at  a  shock  wave  peak  pressure  of  1.00 
kilobar. 


Salinity 
parts  per  1000) 


Temperature 


32 

32 
35 


15 
23 
15 


7.81 
7.31 
7.76 
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by  log(v£0,  S]/?Cp,  ■£])>  an^  -)  the  calculation 
of  quantities  defined  by  line  integrals  along 
aciiabatics  is  greatly  simplified  by  taking  5 
instead  of  T  as  one  of  the  independent  variables. 
The  function  .4 [5]  is  related  simply  to  cu  the 
sound  velocity  at  zero  pressure  and  entropy  5 
according  to  Eq.  (.2.4)  as  follows 

^[S]  =  ci5/"fi;  Ci=*c[Q,  S].  (2.14) 
On  the  basis  of  Bridgman's  p-v-T  data  for  pure 
water,  an  average  value  of.  n  equal  to  7.15  has 
been  selected  for  the  present  calculations.  In 
Section  2  of  Appendix  II,  it  is  shown  that  n 
deviates  from  this  value  by  less  than  4  percent 
in  a  large  pressure-temperature  held  bounded  by 
adiabatics  starting  at  zero  pressure  and  tem- 
peratures of  20°C  and  60°C,  respectively,  and 
extending  to  pressures  of  25,000  kg/cm1.  We 
assume  that  n  has  the  same  value  for  an  aqueous 
0.7  molai  NaCl  solution  as  for  pure  water,  and 
we  obtain  by  interpolation  the  required  values 
of  B(t)  from  R.  E.  Gibson's  values  of  B{t)  for 
dilute  aqueous  NaCl  solutions  (see  Appendix  II, 
Section  1).  The  appropriate  heat  capacity  and 
thermal  expansion  data  are  discussed  in  Section 
1  of  Appendix  II. 

We  now  proceed  to  the  calculation  of  the 
quantities  u,  U,  c,  <r,  and  w.  We  first  express  these 
quantities  with  use  of  Eq.  (2.12)  in  terms  of  p, 

v  =  v{p,  T)  -rfj),  S],     vo=»tf(0,  To), 

t>i-«(0,  TJ=>v[p,  S],  and  d=»c(0,  7*0  -c[0,  5] 
(see  Eqs.  (2.1)-(2.6),  also  Fig.  1)  as  follows: 


M  =-[>(».  -»)]4, 

(2.15) 

U^pva/u, 

(2.16) 

c»Ci(ffi/V)(~l)'2, 

(2.17) 

2c, 

<r= 0./»)(-n'*-l], 

(2.18) 

*=— C(^)"--l], 

(2.19) 

71-1 

Once  the  temperature  7\,  to  which  an  element 
of  fluid  returns  along  the  adiabatic  intersecting 
the  Hugoniot  curve  at  (p,  T),  is  determined,  all 
of  the  above  quantities  may  be  determined  as 
functions  of  p.  To  accomplish  this,  the  enthalpy 
increment,  1H,  occurring  in  the  third  Hugoniot 
condition,  Eq.  (2.3),  is  written  as  the  sum  of  two 
line   integrals,   the  first  along  an   isobar  from 


0.  To)  =  ~_Q.  5,]  to  ;0.  7*0  =  [0,  5]  and 
second  along  an  idiabatic  from  '0,  J'j  =  i~ 
to  tp,  T)=[j>,S2  'see  Fig.  1),  giving: 

oj=  I    vip',  SZdp', 

k=  f    7T0,  S'A4S' -  f     c,(0,  T)iT, 


( 


where  «  is  the  undissipated  enthalpy  air 
defined  by  Eq.  (2.6)  with  po=*Q  and  givei 
plicitly  in  terms  of  ct,  v\,  and  v  in  Eq.  (2 
The  dissipated  enthalpy  h  can  be  determin* 
an  explicit  function  of  To  and  Ti  from  sp< 
heat  data  (Appendix  1 1 ,  Section  1 ) .  Combinin 
third  Hugoniot  condition,  Eq.  (2.3),  with 
(2.19)  and  (2.20)  we  obtain  the  relation 

h      lr      »+i  i 

—  m y (jU-i/»»_l)_v-i/» 

ds     2nL       n  — 1 


V\—  Vo 


2nv\ 

Table   IV.   Properties  of  sea  water  at  a  shock 
(Initial  temperature  09C;  salinity  0.7  m  NaCl;  Co' 
m/sec) 


p 

■ 

U 

c 

3 

«xio-« 

h 

'illo 

(m/ 

(m/ 

Cm/ 

(m/ 

(m/ 

joule/ 

bar) 

sec) 
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MC) 

vex.'- 

SBC)* 

imi 

0 
5 

10 
IS 
20 

:s 

JO 
35 
40 
SO 
60 

:o 
so 

90 


0 
2S7.0 
433.0 
57S.O 
697.5 
SOS  .5 
905.0 
997.0 

1080 

1240 

1385 

1515 

1635 

1740 


1930 
2290 
2585 
2845 
3075 
3235 
3480 
3665 
4000 
4300 
4585 
4855 
5120 


2190 
2720 
3145 
35 10 
3S35 
4125 
+395 
4640 
5095 
5495 
5870 
6225 
6570 


0 
253.5 
+20.5 
552.0 
664.0 
763.5 
355.5 
940.5 

1020 

1175 

1J15 

1+55 

1585 

1705 


0 

0+565 

0.8720 

1.270 

1.655 

2.030 

2.405 

2.770 

3.1+0 

3.860 

+.575 

5.285 

6.000 

6.730 


0 

6.740 
25.80 
54.+0 
36.55 

122.5 

160  J 

201.5 

244.0 

331.0 

419.0 

5090 

595.S 

676.0 


Table  V.   Properties  of  sea   water  at  a   shock 
(Initial  temperature  20°C;  salinicy  0.7  m  NaCl;  C.= 

m/sec.) 


p 

H 

u 

c 

9 
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h 
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(m/ 

(m/ 

(m/ 
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Cm/ 
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bar) 

sec) 

sec) 
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0 
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0 

0 

0 

s 
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io 
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233S 

27S5 
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23.+S 

15 
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3175 
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+9.35 

20 

689.0 
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30.0S 

25 

79J.O 
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3855 

T65.0 
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115.0 

30 

397.5 

3320 

4140 

359.0 

2.+2S 

1S2.S 

35 

990.0 

3510 

+405 

946.5 

2.795 

192.0 

+0 

1075 

3690 

+650 

1030 

3.160 

233.0 

50 

1235 

+020 

5100 

MS5 

1.885 

317.5 

60 

1380 

+325 

5505 

1330 

+.605 

+04.5 

70 

1510 

4610 

5880 

1+65 

5.520 

+89.S 

30 

1625 

4885 

6240 

1600 
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573-S 
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Table  VI.  Properties  oi  ?ea  water  at  a  shock  front. 
Initial  temperature  40=C;  salinity  0.7  m  N'aC!., 


? 

u 

r- 

9 

jXIO"' 

h 

5 

kilo- 

'ml 

fm/ 

•m/ 

m  1 

m. 

joule/ 
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W) 
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?mi 

0 

0 

— 

— 

) 

0 

0 
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10 
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2360 

2775 

+  15  0 
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51 
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31  OS 
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1  290 

+8.30 
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20 

089  0 

2900 

SS50 
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78.70 

.7621 

25 
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3130 
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7+41 

10 
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3335 
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2. +45 

151.0 
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35 
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3325 

+415 

95S.5 

2.315 

139.5 

7179 

40 
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3705 

+<x>0 

1040 
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230.5 

.7080 

50 

1240 

+035 

5110 
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3.915 

315.0 

6926 

00 

1380 

+3+0 

5515 

1345 

+  .040 

+00.5 

0813 

ro 

1510 

+635 

5900 

1485 

5J70 

+85.0 

.6737 

where  v  =  <ui/'u) ".  With  the  aid  oi  tables  of  n/ci* 
and  v\  as  functions  of  7\  and  T0,  Eq.  (2.21)  may 
be  solved  by  successive  approximations  giving  7\ 
as  a  function  of  the  parameter  y.  Since  the  equa- 
tion of  state,  Eq.  (2.12),  may  be  expressed 
simply  as  p  =  B(T:  —  273.16)Qy  —  1],  the  tem- 
perature T\  may  be  determined  as  a  function  of 
the  pressure  p  by  a  tabular  elimination  of  y.  By 
graphical  interpolation,  T-,  is  finally  determined 
for  the  desired  integral  values  of  p  (in  kilobars) , 
and  the  functions  u,  U,  c,  <x,  and  «  are  then  com- 
puted as  functions  of  p  by  means  of  Eqs.  (2.15)- 
(2.19). 

LTL  NUMERICAL  RESULTS  OF  ARONS  AND 
HAL  VERSO N 

In  fundamental  shock  wave  studies,  it  is  fre- 
quently necessary  to  know  values  of  U  —  c0/co 
and  u  at  given  pressure  levels  to  the  highest 
possible  degree  of  accuracy.  With  this  object  in 
view,  the  calculation  methods  described  in 
Section  2A  were  applied  to  the  best  available 
equation-to-state  data.  The  numerical  results  are 
given  in  Tables  I  and  II.  A  critical  discussion  of 
the  equation-of-state  data  will  be  found  in  Ap- 
pendix I  together  with  references  to  the  sources 
from  which  they  were  obtained. 

Table  I  gives  results  for  the  "low  pressure'' 
region,  covering  shock  wave  peak  pressures  of 
from  0  to  1.50  kilobars  {ca.  22,000  p.s.L).  The 
calculations  in  this  table  were  based  upon  the 
Ekman  equation-of-state  for  sea  water  (see  Ap- 
pendix I)  which  is  used  in  the  calculation  of 
sound  velocity  for  echo-ranging  tables. 

Since  the  Ekman  equation  deviates  appreci- 
ably from  experimental  compressibility  data  at 
pressures  exceeding  2  kilobars.  this  equation  was 
abandoned  in  the  "intermediate  pressure"  region. 


The  results  in  Table  II  are  applicable  principalis- 
to  the  region  between  1.5  and  14  kilobars  \ca. 
200.000  p.s.i.)  and  are  based  on  a  careful  fit  of 
the  Tait  equation  to  Adams's  experimental  com- 
pressibility data  (see  Appendix  I). 

Tables  1  and  II  were  computed  for  certain 
specific  values  of  temperature  and  sea  water 
salinity  equivalent  to  0.675  moial  NaCl),  and 
it  is  shown  in  Table  III  that  the  value  of 
U—Cai  Cq  is  not  very  sensitive  to  changes  in  these 
variables. 

IV.  NUMERICAL  RESULTS  OF  KIRK7WOOD 
AND  RICHARDSON 

In  Tables  IV  to  VI,  the  particle  velocity  u,  the 
shock  front  velocity  U,  the  sound  velocity  c,  the 
Riemann  cr-f unction,  the  undissipated  enthalpy 
w,  the  dissipated  enthalpy  h,  and  the  specific 
volume  v  of  sea  water  (0.7  molal  NaCl  solution) 
are  presented  as  functions  of  pressure  p  along 
three  Hugoniot  curves,  starting  at  zero  pressure 
and  the  temperatures  0°C.  20°C  and  40°C, 
respectively.  These  results  have  been  calculated 
by  the  procedures  of  Part  b  of  Section  2  and 
the  data  of  Appendix  II.  The  results  above  30 
kilobars  represent  extrapolations  beyond  the 
range  of  experimental  data;  consequently  the 
validity  of  the  results  above,  say,  50  kilobars,  is 
questionable. 

In  closing  this  discussion  of  the  calculations, 
the  authors  wish  to  acknowledge  their  gratitude 
and  appreciation  to  Professor  J.  G.  Kirkwood  of 
Cornell  University  for  his  contributions  in 
initiating  the  work  and  in  supplying  valuable 
guidance  and  advice. 

APPENDIX  Itt 

1.  Salinity  and  Temperature  Conditions 

All  calculations  were  made  for  sea  water 
having  a  salinity  of  32  parts  per  thousand  (the 
average  salinity  of  sea  water  at  Woods  Hole, 
Massachusetts).  Salinity  is  defined  in  terms  of 
directly  measured  chlorinity  as: 

5  =  0.030  +  1.8050  CI 

where  J  and  CI  are  expressed  in  parts  per 
thousand. 

It  was  calculated  from  the  average  composition 

t"f  Equation-ot-state    data     used     in    computation     oi 

Tabies  I  and  II. 


19 


WATER  AT  THE  FRONT  OF  A  SHOCK  WAVE 


of  sea  water  that  (on  the  basis  of  ionic  strength)  a 
salinity  of  32  parts  per  thousand  is  equivalent 
to  an  NaCl  solution  having  a  molality  of  0.675 
or  a  weight  percentage  of  3.79  percent  NaCl. 

Table  1  was  computed  for  an  initial  tempera- 
ture of  15°C  because  this  temperature  is  a  rough 
average  of  conditions  normally  encountered  in 
experimental  work.  Table  II  was  computed  for 
an  initial  temperature  of  25°C  because  this  was 
the  temperature  quoted  for  the  available  com- 
pressibility data.8  Table  III  shows  that  the 
results  are  not  sensitive  to  small  variations  in 
temperature  and  salinity. 

2.  Specific  Volume  and  Coefficient  and 
Thermal  Expansion 

The  best  sources  of  data  seem  to  be  the 
oceanographical  tables  of  Knudsen.9  Second 
power  equations  in  t(°C)  were  fitted  to  the  data 
tabulated  for  5  =  32: 

For  Table  I : 

v(t)  =0.97709  +  2.05  X10"+*-  IS) 

+4XlO-*(/-15)*. 
For  Table  II: 

v(t)  =0.97956+2.85  X  10~ '(£  -  25) 

+4XlO-«(/-25)-\ 

3.  Heat  Capacity 

The  heat  capacity  data  used  in  computing 
Tables  I  and  II  are  those  quoted  by  S.  Kuwa- 
hara:19 

C,  -  C>  -  0.000422&+0.00000632  Ifi  cal./gmaC, 

Table  VII.  Comparison  of  experimentally  measured 
sound  velocity  with  calculations  based  on  the  Ekman 
Compressibility  Equation.  ;Salinity  — 31.7  parts  per 
thousand.) 


Table  VIII.  Comparison  :>i  \damss  experin 
compressibilities  and  the  empirical  tit  given  by  the  E 
and  Tait  equations. 


Temperature 

Velocity  oi 

»und   •'.  xc 

Deviation 

t°C) 

Measured 

Calculated 

{%) 

10.9 

4887.7 

4875.6 

0.26 

11.6 

+885.8 

4883.8 

0.04 

11.6 

4893.1 

4883.8 

0.19 

11.5 

4902.4 

4882.5 

0.41 

11.1 

4888.3 

4878.5 

0.20 

»  Adams.  J.  Am.  Chera.  Soc.  53,  3769  (1931). 

9  Oceanographical  Tables,  Comissariat  of  Agriculture, 
USSR,  Moscow,  1931.  A  general  compilation  of  oceano- 
graphic  data  by  N.  N.  Zubov.) 

'■'  S.  Kuwahara,  Velocity  of  Sound  in  Sea  Water  and 
Calculation  of  the  Velocity  for  Use  in  Sonic  Sounding 
Hydrographic  Dept.  I.J.N.  Tokyo,  1938). 


?«—»),'»» 

Adams 

tkman 

Tait*  equaci 

experimental) 

equation 

n  -7.300 

p 

Pure 

3.79% 

5-32 

3-3.012        i  ' 

ilc  bar) 

HiO 

NaCl 

.Table  [) 

Table  II)      S- 

0.0 

0.0000 

0.0000 

0.0000 

0.0000         0. 

0.5 

.0212 

.0196 

.0198 

.0197 

1.0 

.0393 

0368 

.0370 

.0368 

1.5 

.0535 

.0522 

.0522 

.0518 

2.0 

.0699 

.0658 

.0655 

.0653 

3.0 

0945 

.0894 

.0871 

.0887 

4.0 

.1152 

.1091 

— 

.1083 

5.0 

.1330 

.1265 

— 

.1254 

6.0 

.1485 

.1417 

— 

.1405 

7.0 

.1622 

.1552 

— 

.1540 

8.0 

.1746 

.1670 

— 

.1662 

9.0 

.1858 

.1781 

— 

.1775 

10.0 

.1964 

.1886 

— 

.1876 

11.0 

.2059 

.1980 

~~ — 

.1972 

• :«»-») /ti«<  i'<i  io«(i  —p/B). 


rhere 


C„-  1.005  -0.004136s+9.0001098s2 

-0.0000013 

In  the  above  equations,  t  is  temperature  i: 
and  s  is  salinity  in  parts  per  1000.  These 
are  in  good  agreement  with  those  used  by  I 
wood  and  Richardson,  quote  in  Appendix  I 

4.  Compressibility  Data  for  Low  Pressure  Re 
(Table  I) 

The  following  equation  was  used  in  compi 
Table  I : 

4886 

10*u  = [227  +  28.33/  -0.551*2 

1+0.183? 

+0.004<,]+p[105.5+9.50*-0.158J 

/r-28x 
-  \.SpH~  I J[147.3  - 2.72/ 

+0.04<2-p(32.4-0.87^0.02*2)] 

/7-28\: 
+  ( J  [4.5-0.1/-p(1.8-0.0 

where  p  is  pressure  in  kilobars,  t  is  tempers 
in  degrees  centigrade,  and  n  is  defined  by: 

:>=t>o(l-up), 
y  is  defined  by : 

7=  -0.069  +  1.4708  CI -0.001570  CI2 

+0.000039* 
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The  above  empirical  equation  tor  sea  '.vater 
compressibility  is  due  to  Ekmana  and  has  been 
widely  used  for  computation  of  sound  velocity 
in  sea  water.1012  The  validity  of  the  Ekman 
equation  for  sound  velocity  calculations  was 
verified  experimentally  as  indicated  in  Table  VII. 

Experimental  sound  velocity  measurements 
were  made  by  recording  with  a  rotating  drum 
camera  the  signals  applied  to  a  cathode  ray 
oscilloscope  by  two  very  small  piezoelectric 
gauges  placed  a  known  distance  apart.  The  sound 
source  was  a  No.  S  detonator  cap  placed  far 
enough  away  from  the  gauges  so  that  the  effect 
of  finite  pressure  amplitude  was  less  than  0.03 
percent.  An  error  of  about  0.2  percent  was  in- 
herent in  the  experimental  work  owing  to  slight 
errors  in  the  alignment  of  the  two  recording 
gauges  with  the  sound  source.  This  accounts  for 
the  magnitude  and  systematic  nature  of  the  dis- 
crepancy apparent  in  Table  VII. 

Further  verification  of  the  applicability  of 
Ekman's  equation  in  the  region  up  to  1.50 
kilobars  is  given  in  Table  VIII  where  values  ob- 
tained from  the  equation  are  compared  with  the 
experimental  values  of  Adams  for  N'aCl  solutions. 

5.  Compressibility  Data  for  Intermediate  Pres- 
sure Region   Table  II) 

As  indicated  in  Part  1  of  this  appendix,  a  sea 
water  salinity  of  32  parts  per  thousand  cor- 
responds to  a  3.79  weight  percent  solution  of 
NaCl.  The  compressibility  of  NaCl  solution  of 
this  concentration  was  obtained  by  graphical 
interpolation  of  Adams'3  data.' 

The  Tait  equation  in  the  form : 

o(0,  T)-v(p,  T)/v(Q,  Z")-U/>*)  log£l+p/B(t)2, 

i  =  (r-273.16)°C 

was  then  fitted  to  Adams's  data.  In  an  effort  to 


Table 

IX.   Val 

fusing 

jes    if 

n»  7.15 

i  computed   from   p  —  v—T  data 
in  computation  of  XT). 

CO 

5000 

p  (kg/cm1) 
15.000 

25.000 

20 
40 
60 

7.211 
7.360 
7.411 

7.183 

7.126 
7.054 

7.130 
6.969 
6.86S 

''••centigrade   temperature   through   wtiich    the  adtabatic  for  5 
passes  at  zero  pressure. 


11 V.  W.  Ekman.  Publications  de  Circonstance  So.  43 
'Conceil  Permanent  Internationale  Pour  L' Exploration  de 
la  Mer.  November  1908). 

11  Matthews,  Tables  of  the  Velocity  of  Sound  in  Pure 
Water  and  Sea  Water  for  Use  m  Echo  Sounding  and  Sound 
Ranging  'Hydrographic  Dept.,  Admiralty,  H.D.  No.  232). 


check  the  Tait  equation  against  the  Ekman 
equation  used  for  computation  of  Table  i.  a  fit 
was  first  made  to  the  lower  pressure  region. 
Values  of  n  and  5(  25CC)  were  so  selected  that  the 
equation  not  only  fitted  the  data  of  Adams  with 
adequate  precision  but  also  yielded  the  correct 
velocity  of  sound  in  the  limit  of  zero  pressure. 
This  additional  restriction  that  the  equation 
give  co=  1528  m  sec.  at  25°C  and  s  =  32)  required 
that  nB(2S°C)  =  23.497,  the  latter  relation  being 
obtained  from  the  thermodynamic  equations: 

(dv  \  iv      T  /  dv  \  ; 

dp/  t        c<r     cp\  dT/  p 


/dv\  va 

—  )    = atp< 

\dp/  r         nB{i) 


-0. 


In  this  case  n  was  taken  as  7.445  -and  3(25°C) 
as  3.156  kilobars,  and  the  resulting  equation  fits 
the  data  of  Adams  quite  closely  up  to  pressures 
of  about  4  kilobars  as  shown  in  Table  VIII.  For 
purposes  of  further  calculation,  the  temperature 
variation  of  B  was  assumed  to  be  the  same  as  that 
used  by  Kirkwood  and  Richardson  on  the  basis 
oi  a  private  communication  from  Gibson  ;see 
Appendix  II).  Calculation  of  U  —  Cg/Co  at  1.00 
kilobar  yielded  a  value  of  7.85  percent,  in  good 
agreement  with  the  value  of  7.81  percent  ob- 
tained from  the  Ekman  equation. 

Having  verified  the  accuracy  of  results  ob- 
tained from  the  Tait  equation  when  fitted  as 
described  above,  the  same  technique  was  used  to 
fit  the  equation  to  the  intermediate  pressure 
range  (up  to -values  for  11  kilobars  quoted  by 
Adams;  it  was  assumed  safe  to  extrapolate  the 
resulting  equation  to  pressures  of  14  or  15  kilo- 
bars). It  was  found  that  the  best  fit  of  the  data 
as  well  as  a  correct  value  for  the  velocity  of 
sound  were  obtained  by  taking  rt  =  7.300  and 
3(25°C)  =3.012,  the  temperature  variatior  of  3 
again  being  assumed  to  be  that  mentioned  above. 
The  Tait  equation  containing  these  parameters 
was  then  used  for  the  computation  of  Table  II. 
The  fit  of  the  equation  to  Adams's  data  is  shown 
in  Table  VIII. 

APPENDIX  H 

1.  Data  Employed  in  the  Computations  of  Part 
B  of  Section  II 

In  the  modified  Tait  equation,  Eq.  (2.12),  the 
function  4 [S]  =»£(*) •  where  i  =  T[0,  5] -273. 16, 
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is  determined  from  the  empirical  values  of  Bit), 
fitting  the  original  isothermal  Tait  equation,  Eq. 
1,2.7),  to  experimental  data.  R.  E.  Gibson11  gives 
third-degree  ^-expansions  of  3d)  for  various 
molalities  of  NaCl.  By  interpolating  the  coef- 
ficients (the  constant  term  numerically  and  the 
other  graphically)  for  a  molality  of  0.7,  one 
obtains  Bit)  =3.134-  1.65  X  10~3(i- 55)  -1.181 
X  lQr*(t-SS)*+S.32 X  10-7^ -55)'  kilobars. 

The  specific  heat  cp(,0,  T)  for  a  0.7  molal  NaCl 
solution  was  obtained  by  interpolation  from  the 
values  quoted  in  the  International  Critical  Tables 
and  Physikalischchemtsche  Tabellen.  The  resulting 
set  of  values  is  fitted  adequately  by  the  ex- 
pression 

c,(0.  t  +  273. 16)  =3.9644+6.24 

X  I0~*t  joule,  gm.  deg. 

From  Gibson  and  Loefner14  a  set  of  values  of 
t/(0,  T)  covering  the  range  from  25°C  to  95°C 
inclusive  was  obtained  for  a  0.7  molal  NaCl 
solution  by  means  of  empirical  equations  giving 
v(0,  T)  as  a  function  of  concentration  for  each 
temperature.  In  extrapolating  to  higher  tem- 
peratures, the  relation, 

»(0,  *  +  273.16)=0.994150+2.929Xl0-'(;-25) 

+  3.241  X  10-*((-25)2  cm'/gm. 

was  used;  for  lower  temperatures  (f<10°C), 

v(0,  f+273.16)  -0.991442+6.025 

Xl0-«(*-3.8)2cmJ/gm. 

2.  Test   of   the   Modified   Tait   Equation   with 
Bridgman's  Data  for  Pure  Water.  Deter- 
mination of  the  Characteristic 
Constant  n. 

The  modified  Tait  equation-of-state,  Eq. 
(2.12),  may  for  our  present  purposes  be  written 
in  the  form 

log(t>[0,  5]/vfj>,  SJ) 

=-(l/n)log(l+pM[5])     (II-l) 

where  A £S]  is  related  to  the  Bit)  in  the  original 
isothermal  equation  of  state  as  follows, 

<•-r.-273.16.  (II-2) 

ra=r[o,s]. 

"  Private  communication. 

14  Gibson  and  Loerfler,  J.  Am.  Chem.  Soc.  S3,  443  0941). 
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According  to  the  convention  introduced  in 
B  of  Section  II,  parentheses  (  )  after  a  fur 
denotes  that  the  independent  variables  are 
r,  whereas  square  brackets  [  ]  denote  that 
are  p  and  5. 

Now  we  wish  to  test  Eq.  (II-l)  with  I 
man's15  p-v-T  data  for  pure  water  witl 
ultimate  object  of  finding  the  best  value 
We  assume  implicitly  that  n  does  not 
rapidly  with  NaCl  concentration.  To  mak 
comparison,  we  first  must  know  the  values  i 
temperature  T  corresponding  to  the  vj 
points  [_p,  S~\,  the  calculation  of  which  we 
sider  below. 

Letting  HO,  5]  =  r«,  Tfp,  5]  -  T,  and  2 
=  Ar,  we  have 


Ar  = 


I 


'  <37Tj>.  S]  r"  *»L>.  S] 

dp  dS 


-dp 


-f 


■dp-     ( 


Using  Eq.  (II-l),  a  simple  calculation  yielc 
G 


*T> Z(l-rDKl+p'A) 

(l+p/A)l'n 

-(»-Z?)(l+^//l)1"'+n-l]>     ( 


where 


A»A£S1-B(t0), 

A'ZSy[0,S]     roB'(h)-v(Q,To) 


nA£S} 


£=• 


n-1  (»-l)-c„(0.  To) 

dv(Q,  5] 

dS  nBit0)do 


/I'CSXO.  5]      B'MviO,  To) 
'dv{0,  21' 


/dv{0,  f)\ 

V      dT     /t-To 


To  calculate  T,  given  a  specified  p 
r0—  rfO,  S],  a  tentative  value  of  »  —  7.15 
chosen  for  use  in  Eq.  (1 1-4).  The  correspor 
value  of  v£p,  S]  =  vip,  T)  was  obtained  by  i 
poiation  from  Bridgman's14  p-^o-T  data, 
serting  these  values  of  v\j>,  5]  in  Eq.  (II-l) 
knowing  the  values  of  v[0,  5]=v(,0,  To) 
Bit0)  for  pure  water,  a  set  of  values  of  n 
calculated  for  p  =  5,000,   15,000,  25,000  kg 

"Bridgman.  J.  Chem.  Phys.  3,  597  (,1935)  and  p 
communication. 
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and  ;0=  To -273. 16  =  20°,  40°,  60°C.  the  results 
summarized  in  Table  IX. 

Additional  data  Jor  pure  water;  used  in  Eqs. 
TI-1)  and  (II-4)  were 

3(i)  =  2.996+7.285  X  10~3'i  —  25)  -  1.790 

XlQ-*(t-2S)i-r6A3Xl0-:(t-2S)1    [ritobars,1* 

and 
1      d  lognO,  273.16+0 


.303 


St 


2(*-3.98) 


244.860+  15,040  (J-  3.98)°" 

(0.62)(15,040)U-3.98)'« 
[244,860+ 15, 040(/  -  3.98)0*2]2 

obtained  from  Ipatov's14  empirical  equation  for 
v  by  differentiation. 

The  average  value  of  n  is  7.146.  In  the  present 
calculations  this  value  has  been  rounded  off  to 
7.15. 

The  entries  in  Tables  IV,  V,  and  VI  therefore 
contain  more  significant  figures  than  the  test 
justifies..  On  the  basis  of  the  test,  the  errors 
associated  with  the  use  of  the  modified  Tait 
equation  are  of  the  order  of  several  percent.  In 

'•  I.  V.  Ipatov,  J.  Phys.  Chem.  (U.S.S.R.)  5.  1230  (1934). 


particular,  the  results  obtained  for  low  pressures 
will  disagree  with  known  data  by  several  percent. 

APPENDIX  m 

Symbols 

.4 [5]  =  parameter  in  modified  (."adiabatic')  Tait  equauon- 
of-state. 
5(0  =  parameter  in  isothermal  Tait  equation-of-state. 
c  » local  velocity  of  sound. 
Co  =*  velocity  of  sound  at  zero  pressure. 
<;„=■  specific  heat  at  constant  pressure. 

c„(,0.  T'  dT'< 


n  * 
Pi* 

?■■ 
S* 

T-- 
ii  = 


■enthalpy  increment:  \H  «»«+&. 

'characteristic  constant  in  Tait  equation-of-state. 

■initial  pressure  ahead  oi  shock  front,  £o  =  0,  in 
these  calculations. 

pressure  behind  shock  front. 

entropy. 

sea  water  salinity. 

'temperature  in  'C. 

absolute  temperature. 

particle  velocity  behind  shock  front. 

shock  front  propagation  velocity. 

specific  volume  of  medium  ahead  of  shock  tront. 

-specific  volume  oi  medium  behind  shock  front. 

mean  compressibility  at  zero  pressure  over  temper- 
ature range  A7*. 

density. 


a  "  Riemann  function 


'»  cfj/.  5] 

JV     _ 
sip',  Sjip'. 
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APPENDIX  3 

COMPUTER 

LISTING 

01  LBLTHH 

26 

STO  05 

51 

X 

76 

- 

101 

X<Y? 

02  TLOOP? 

27 

TPRESSURE? 

52 

RCL  25 

77 

STO  07 

102 

GTCTSH0C] 

03  PROMPT 

28 

PROMPT 

53 

- 

78 

T 
HH  =» 

103 

GTO  01 

04  STO  34 

29 

STO  12 

54 

CHS 

79 

ARCL  X 

104 

LBl/VW 

05  TN? 

30 

RCL  05 

55 

RCL  25 

30 

AVIEW 

105 

25 

06  PROMPT 

31 

XEQTVW 

56 

RCL  23 

81 

STOP 

106 

- 

07  STO  28 

32 

STO  27 

57 

YtX 

82 

RCL  06 

107 

STO  03 

08  1 

33 

LBL  01 

58 

- 

33 

1000 

108 

SNTER+ 

09  + 

34 

XEQ1! 

59 

2 

34 

X 

109 

X 

10  RCL  28 

35 

STO  25 

60 

/ 

35 

T 
HT  =» 

110 

3.241E-6 

11  1 

36 

RCL  01 

61 

RCL  20 

36 

ARCL  X 

111 

X 

12  - 

37 

XEQTVW 

62 

X 

37 

AVIEW 

112 

RCL  03 

13  / 

38 

STO  26 

63 

RCL  26 

88 

STOP 

113 

2.929E-4 

14  STO  21 

39 

RCL  01 

64 

RCL  27 

39 

RCL  07 

114 

X 

15  1 

40 

XEQTSPEED 

65 

- 

90 

RCL  06 

115 

+ 

16  RCL  28 

41 

STO  20 

66 

2 

91 

1000 

116 

0.99415 

17  1/X 

42 

RCL  01 

67 

/ 

92 

X 

117 

+ 

18  - 

43 

XEQTHT 

63 

RCL  20 

93 

- 

118 

RTN 

19  STO  22 

44 

STO  06 

69 

X 

94 

THH-HT  = 

119 

L3LTY 

20  RCL  28 

45 

RCL  25 

70 

RCL  26 

95 

ARCL  X 

120 

TTEMP  1? 

21  1/X 

46 

RCL  22 

71 

/ 

96 

AVIEW 

121 

PROMPT 

22  CHS 

47 

Y+X 

72 

RCL  25 

97 

STOP 

122 

STO  01 

23  STO  23 

48 

1 

73 

1 

98 

ABS 

123 

XEQT3T 

24  TTEMP  0? 

49 

- 

74 

- 

99 

RCL  34 

124 

RCL  08 

25  PROMPT 

50 

RCL  21 

75 

X 

100 

X<>Y 

125 

1/X 

25 


126  RCL  12 

127  X 

128  1 

129  + 

130  STO  25 

131  RTN 

132  LBLTS?EED 

133  STO  01 

134  XEQTBT 

135  STO  11 

136  RCL  01 

137  XEQTVW 

138  0.001 

139  X 

140  RCL  11 

141  X 

142  STO  20 

143  RTN 

144  LBLTHT 

145  STO  01 

146  RCL  05 

147  - 

143  3.9644 

149  X 

150  RCL  01 


151  ENTERt 

152  X 

153  RCL  05 

154  ENTERt 

155  X 

156  - 

157  0.000312 

158  X 

159  + 

160  RTN 

161  L3LTBTPRIME 

162  ENTERt 

163  55 

164  - 

165  STO  10 

166  -2.362E-4 

167  X 

168  -0.00165 

169  + 

170  RCL  10 

171  ENTERt 

172  X 

173  1.596E-6 

174  X 

175  + 


176  1.013E8 

177  X 

173  STO  24 

179  RTN 

180  L3LTBETA 

181  25 

182  - 

183  6.482E-6 

184  X 

185  2.929E-4 

186  + 

187  0.001 

188  X 

189  STO  19 

190  RTN 

191  L5LTBT 

192  55 

193  - 

194  ENTERt 

195  ENTERt 

196  STO  04 

197  X 

198  STO  02 

199  X 

200  5.32E-7 


201  X 

202  RCL  02 

203  -1.131E-4 

204  X 

205  + 

206  RCL  04 

207  -0.00165 

208  X 

209  + 

210  3.134 

211  + 

212  STO  08 

213  1.013E8 

214  X 

215  STO  09 

216  RTN 

217  LBLTSH0CK 

218  RCL  11 

219  1.013E8 

220  / 

221  1/X 

222  RCL  12 

223  X 

224  1 

225  + 


26 


226  RCL  28 

227  1/X 

228  Y+X 

229  1/X 

230  RCL  26 

231  X 

232  STO  13 

233  FIX  6 

234  TV  =» 

235  ARCL  X 

236  AVIEW 

237  STOP 

238  FIX  0 

239  RCL  12 

240  1.013E8 

241  X 

242  STO  15 

243  RCL  27 

244  RCL  13 

245  - 

246  0.001 

247  X 

248  X 

249  SQRT 

250  STO  14 


251  "WATER  0 

252  ARCL  X 

253  AVIEW 

254  STOP 

255  RCL  15 

256  RCL  27 

257  0.001 

258  X 

259  X 

260  RCL  14 

261  / 

262  STO  16 

263  T5H0CK  U 

264  ARCL  X 

265  AVIEW 

266  STOP 

267  RCL  20 

268  RCL  28 

269  X 

270  SQRT 

271  STO  17 
27  2  RCL  26 

273  RCL  13 

274  / 

275  RCL  28 


276  1 

277  - 
273  2 

279  / 

280  Y+X 

281  X 

282  STO  18 

283  WAVE  C  = 

284  ARCL  X 

285  AVIEW 

286  STOP 

287  GT0TDEL  T 

288  RTN 

289  LBLTCP 

290  0.624 

291  X 

292  3964.4 

293  + 

294  STO  29 

295  RTN 

296  LBLTG1 

297  RCL  05 

298  273.16 

299  + 

300  RCL  24 


301  X 

302  RCL  27 

303  0.001 
204  X 

305  X 

306  RCL  28 

307  1 

308  - 

309  / 

310  RCL  29 

311  / 

312  STO  30 

313  RTN 

314  L3LTD2 

315  RCL  28 

316  RCL  08 

317  X 

313  RCL  19 

319  X 

320  RCL  24 

321  1.013E8 

322  / 

323  / 

324  RCL  27 

325  0.001 
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326  X  351  RCL  31 

327  /  352  + 

328  STO  31  353  X 

329  RTN  354  CHS 

330  LBLTDEL  T  355  RCL  28 

331  RCL  12  356  + 

332  RCL  08  357  1 

333  /  358  - 

334  1  359  RCL  32 

335  +  360  RCL  31 

336  STO  32  361  1 

337  RCL  05  362  + 

338  XEQTCP  363  X 

339  RCL  05  364  + 

340  XEQTBTPRIME    365  RCL  30 

341  RCL  05  366  X 

342  XEQTBETA  367  RCL  33 

343  XEQTG1  368  / 

344  XEQTD2  369  FIX  2 

345  RCL  32  370  TDEL  T 

346  RCL  28  371  ARCL  X 

347  1/X  372  AVIEW 

348  Y+X  373  RTN 

349  STO  33  374  END 

350  RCL  28 
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APPENDIX  C 

CALCULATION  OF   TEMPERATURE,    INTERNAL   ENERGY,    AND   ENTHALPY 
FOR   ISENTROPIC  FLOW 


In  order   Co  calculate   the  flow  field  downstream  of   the  shock  wave, 

as   illustrated   in  Figure  C-l,    one  needs    to  calculate  the   internal   energy 

and   enthalpy  of    the  water   as   it   expands.      The  energy  equation  appropriate 
for    the   flow  illustrated   in  Figure  C-l   is 


3L                         i.  r,                                                             < 

r          /               U        +    V     .  ,             d      r          ,               U 

3~{pu(e  + = )  +  pu]  +  f-[pv(e  +  - 

ox                         -  or 


")   +  pv] 


2         2 
+  ^[pv(e  +  u     \  V   )   +  pv]    -  0 


(C-l) 


Equation   (C-l) ,   which  assumes   inviscid   flow  without  heat  conduction,    indicates 
the  need   to  calculate  internal   energy. 


SHOCK  WAVE 

NUMBERS  TO 
INDICATE  LOCATION 
ALONG 
STREAMTUBE 


ALUMINUM  JET 


Figure  C-l.   Sketch  of  Flow  Field  Due  to  an  Aluminum  Jet  Penetrating  Water, 
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A.   Calculation  of  Internal  Energy  and  Enthalpy 

3y  definition  internal  energy  and  enthalpy  are  related  by 

H  -  e  +  pv  (C-2) 

Capital  H  is  used  for  enthalpy  since  h  is  used  for  dissipated  enthalpy. 
Refer  to  Figure  C-2  which  is  similar  to  Figure  1  of  Appendix  A.   The 
points  indicated  in  Figure  C-l  are  shown  also  in  Figure  C-2.   For 
example,  in  Figure  C-l  point  3  is  immediately  downstream  of  the  shock  wave. 
In  passing,  one  should  note  that  to  calculate  shock  properties,  the 
normal  component  of  velocity  should  be  used.   Point  3  is  shown  in  Figure  C-2, 
As  water  flows  along  the  streamtube,  expansion  occurs.   The  pressure 
decreases  as  indicated  by  point  4  in  Figure  C-l  or  Figure  C-2. 


temperature,  T 


^ 


CONSTANT  PRESSURE 


CONSTANT  TEMPERATURE 
3 


ADIABATIC  (ISENTROPIC) 
KUGONIOT 


pressure,  p 


Figure  C-2.     Hugoniot  and  Adiabatic    (Isentropic)    Curves   in 
the  Pressure-Temperature  Plane. 
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From  equation  2.20  of  Appendix  A 

~H  "  H3  "  H0  "  U31  +  h10  ((>3) 

The  subscripts  31,  used  above  and  later  in  other  equations,  mean  u> 
is  evaluated  along  the  adiabatic  curve  between  points  3  and  1. 
Likewise  h-»  is  evaluated  along  a  line  of  constant  pressure  between 
points  1  and  0  in  Figure  C-2.   From  equation  (2.19)  of  Appendix  A 

c2    v  n~1 
"31-n~T^)    -«  (C"4) 

Also,  from  equation  (2.20)  of  Appendix  A, 

1 

hin  -  /  C  (0,T')dT'  (C-5) 

10   0  P 

At  some  intermediate  point,  such  as  point  4  in  Figure  C-2,  one  replaces 
subscript  3  by  4  in  equations  (C-4)  and  (C-5) . 

Combining  equations  (C-2)  through  (C-5),  an  expression  for  internal 
energy  e  can  be  derived;  the  result  is 

?    n-1  , 

e  =  h.n  +  cA-r TT  +  ~ ^  (c~6) 

10    1  n(n  -  1)    nz   n  -  1 

where  z  is  v  /v.   Hence  equation  (C-6)  uses  v  as  the  independent  variable. 

An  alternate  form  can  be  derived  using  p  as  the  independent  variable;  the 

result  is 

n-1 

2r(l  +  x)  n   .     1         J.,    cr   7, 

e  »  h.n   +  c.  [ — t -r-r—     +  —, -J    (C-7) 

10    1  n(n  -  1)     n(1  +  x)l/n   n-1 

where  x  ■  p/B  ■  p/A.   The  symbols  3  and  A  are  defined  on  page  22,  Appendix  A; 
see  equation  (II-2).   Equations  (C-6)  and  (C-7)  have  a  very  similar  form. 
An  universal  equation  can  be  written  as 

?  n~"^"     i      i 

e  -  h.n  +  eft-— rr-  + ~r]  (C-8) 

10    1  n(n  -  1)   nw   n-1 
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where 

z  for  v  /v  as  independent  variable 
w  =  {       j,  (C-9) 

(1  +  x)    for  p/B  as  independent  variable 

Note  that  h. n  is  a  constant  value  for  any  streamtube. 

3.   Calculation  of  Temperature 

If  the  energy  equation  incorporates  thermal  conduction,  a  term 

involving  temperature  is  used.   The  term  has  the  form 

Q  -  k~  (C-10) 

ox 

2 
where  Q  is   the  heat  flux  in  Joule/m     sec,   and   k  is   thermal  conductivity. 

To  calculate  the   temperature  at  some  intermediate  point,    i.e.   T,, 
one  can  use   the  subroutine  DEL  T.      In  DEL  T,    p-   is  replaced  by  p, .      As 
a  note  of  caution,   when  intermediate  values  of   temperature  are  calculated, 
the  pressure  behind   the  shock  wave,   p_,   is  replaced   in  storage  register  12 
by  p.,    the  intermediate  pressure. 
C.      Subroutine  E-BAR  and  User's  Guide  to  E-BAR 

A  subroutine  has  been  written  to  calculate  internal  energy  and 
enthalpy.     A  listing  of   the  subroutine  is  given  below: 


01  LBLTE-BAR  11  YtX  21  - 

02  TISEN  P?  12  STO  35  22  / 

03  PROMPT  13  RCL  28  23  RCL  28 

04  STO  36  14  1  24  1/X 

05  RCL  08  15  -  25  RCL  35 

06  /  16  YtX  26  / 

07  1  17  RCL  28  27  + 

08  +  18  /  23  RCL  28 

09  RCL  28  19  RCL  28  29  1 

10  1/X  20  1  30  - 
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31  1/X 

41 

AVIEW 

51 

RCL  35 

32  - 

42 

STO  37 

52 

/ 

33  RCL  17 

43 

STOP 

53 

RCL  37 

34  ENTER" 

44 

RCL  36 

5u 

- 

35  X 

45 

1.013E8 

55 

*H  BAR 

36  X 

46 

X 

56 

ARCL  X 

37  RCL  07 

47 

RCL  26 

57 

AVIEW 

38  + 

48 

0.001 

58 

STO  38 

39  TE-BAR  - 

49 

X 

59 

END 

40  ARCL  X 

50 

X 

As  a  convenience,    subroutine  E-BAR  should   be  assigned   to   some  key  when 
calculator  is   in  USER  mode. 

To  operate  E-BAR,    one  must  execute  HH  as  discussed   on  pages  4   to   9 . 
Operation  of  HH  fills   the  registers  with  values  which  are  recalled  for 
use  by  E-BAR.      For  a  sample  problem,   HH  was   executed  using  n  ■  7.15, 
p  ■  90  kilobars,   and   t^  ■  0.     When  HH  has  been  completed,    one  proceeds  as 
follows : 

1.  Continue  with  HP41CV  in  USER  mode. 

2.  Assign  E-BAR  to  a  key  on  the  keyboard  by  pressing 

|~I      ALPHA     E-BAR     ALPHA,    R+ 
In  this   example,    E-BAR  has  been  assigned    to   the  R4-  key. 

3.  Press  R+,   and  observe  ISEN  P?.      This   is   the  isentropic  pressure 
at  point  4  in  Figure  C-2;    obviously  one  can  select  any  value. 
The  number   in  parentheses   is   the  value  used   in   the  sample 
problem   (60  kilobars) . 

4.  Press  R/S. 
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5.  Calculator  displays  E-BAR  3.   For  the  sample  problem 
(E-BAR  =  1,249,766  Joule/kg). 

6.  Press  R/S. 

7.  Calculator  displays  H  BAR  ■.   For  the  example  problem 
(H  BAR  -  5,449,312  Joule/kg). 

Subroutine  E-BAR  has  now  been  completed.   If  calculations  are  desired  at 
other  values  for  p, ,  continue  as  follows: 

8.  Press  R/S,  and  observe  ISEN  P? . 

9.  Input  ISEN  P. 

10.  Press  R/S. 

11.  Etc. 

D..  Subroutine  ISEN  T  and  User's  Guide  to  ISEN  T 

A  subroutine  ISEN  T  has  been  written  to  calculate  temperature  along 
the  adiabatic  (isentropic)  curve  of  Figure  C-2.  A  listing  of  the  subroutine 
is  given  below: 

01  LBLTISEN  T  06  STO  39  11  TISEN  T  - 

02  TISEN  P?  07  RCL  01  12  ARCL  X 

03  PROMPT  08  +  13  AVIEW 

04  STO  12  09  273.16  14  END 

05  XEQTDEL  T  10  + 

To  operate  ISEN  T,  one  must  execute  HH  as  discussed  on  pages  4-9 .  For 
a  sample  problem,  HE  has  been  executed  with  n  a  7.15,  p  ■  90  kilobars,  and 


co-°- 


1.  Continue  with  HP41CV  in  USER  mode. 

2.  Assign  ISEN  T  to  a  key  on  the  keyboard  by  pressing 

|~j   ALPHA  ISEN  T  ALPHA,  EEX 
In  this  example,  ISEN  T  has  been  assigned  to  EEX  key. 
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3.  Press  EEX,  and  observe  ISEN  ??.   This  is  the  isentropic 
pressure  at  point  4  in  Figure  C-2.   Since  o.  or  ISEN  P  is  an 
independent  variable,  any  value  can  be  selected.   The  number  in 
parentheses  is  the  number  used  in  the  sample  problem  (60  kilobars; . 

4.  Press  R/S. 

5.  Calculator  displays  ISEN  T  =».   For  the  sample  problem,  the  result 
is  (ISEN  T  •  740.6).   The  dimensions  are  °K. 

Subroutine  ISEN  T  has  now  been  completed.   If  calculations  are  desired 
at  other  values  for  p.,  continue  as  follows: 

6.  Press  R/S,  and  observe  ISEN  P?. 

7.  Input  ISEN  P  (50  kilobars). 

8.  Press  R/S. 

9.  Observe  ISEN  T  »  691. 7°K. 
10.   Etc. 

E.   Additional  Assignment  of  Storage  Registers  for  Subroutines  E-BAR  and  ISEN  T 
Table  C-I  gives  the  additional  assignment  of  storage  registers. 


Table  C-I. 

Additional 

Storage  Registers  for  E-BAR  and  ISEN  T 

Register 

Symbol 

Definition         Units   Programs 

35 

36 

37 
38 


(1  +  p/B) 


1/n 


e 
H 


Also  symbol  w  in 
equation  (C-8) . 
Equals  v  /v. 

Intermediate 
pressure  on 
isentropic  line. 

Internal  energy. 

Enthalpy. 


39  AT  (Isentropic)  Change  in  tempera- 
ture along  the 
adiabatic  line  in 
Figure  C-2. 


E-3AR 


kilobars   E-BAR 


Joule/kg   E-3AR. 

Joule/kg   E-3AR 
o. 


:< 


ISEN  T 
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F.   Sample  Results  for  Properties  Along  Adiabatic  (Isentropic)  Line 

Table  C-II  summarizes  the  results  for  conditions  along  the  adiabatic 
line  in  Figure  C-II.   Input  values  for  the  shock  wave  were  n  a  7.15,  p  =  90 
kilobars,  t  =»  0°C,  and  t  =»  167.762°C. 


Table  C-II. 

Sample  Results  for 
Adiabatic  (Isentrop 

Conditions 
ic)  Line 

Along 

P 
kilobars 

e 
kilo joule/kg 

H 

kilojoule/kg 

T 
°K 

0 

674 

674 

441 

5 

707 

1180 

460 

10 

755 

1634 

486 

20 

357 

2474 

539 

30 

958 

3260 

591 

40 

1057 

4013 

642 

50 

1154 

4741 

692 

60 

1249 

5449 

741 

70 

1343 

6142 

789 

30 

1435 

6820 

836 

90 

1526 

7488 

883 

For  p  =  0  kilobars,  the  temperature  T  is  the  residual  due  to  entropy  increase 
across  the  shock  wave.   For  p  ■  90  kilobars,  the  temperature  T  is  the  value 
immediately  behind  the  shock  wave. 
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APPENDIX  D 
(NPS67-82-001) 

CALCULATION  OF  PRESSURE  AS  A  FUNCTION  OF  SPECIFIC  VOLUME  AND  INTERNAL  ENERGY 

Introduction 

In  order  to  use  efficiently  the  computer  program  developed  by  Kutler  and 
Chakravarthy  [D-l],  the  ability  to  calculate  pressure,  p,  as  a  function  of 
specific  volume,  v,  and  internal  energy,  e,  is  needed.   The  program  for 
HP41CV,  which  is  described  in  this  appendix,  calculates  p  =  p(v,e). 
Relevant  Equations 

The  relevant  equations  for  the  calculation  can  be  assembled  from  the  equations 
in  this  report.  Define 

z  =  v  /v  (D-l) 


wher 


e  v.  is  the  specific  volume  at  t.. ,  and  v  is  the  specific  volume  at  the 


desired  condition.   The  value  for  v..  can  be  calculated  from 

v(t)  =  0.99415  +  0.0002929(t  -  25)  +  3.241  x  10_6(t  -  25)2  (D-2) 

The  internal  energy  is  given  by 


e  =  hin  +  eh-/ TY  +  i7-^-T]  =  e(t.  ,v  /v)  (D-3) 

10    1  n(n  -  1)   nz   n  -  1       11 

Note  that  e  is  a  function  of  t  and  v. /v.  The  quantity  h^Q   is 


h       =     f  "   c  (0,t)dt  =  3.9644(t_  -  tn)  +  3.12  x  10  4(t?  -  t2) 
10         p  i    u  j.    u 

0 


(D-4) 


The  temperature,    t„,    is   the  reference  temperature  which  can  be  the   initial 
water   temperature  before  the  shock  wave  arrives. 

The  pressure  can  be  related   to  specific  volume  by  the  Tait  equation 

p  =  BCt^MCVj/v)11  -  1]    =  pCv,^)  (D-5) 
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(D-6) 


The  parameter  B(t)  is  given  by 

B(t)  =  3.134  -  1.65  x  10~3(t  -  55)  -  1.181  x  10_4(t  -  55)2 

+  5.32  x  10~7(t  -  55)3 

2 
Finally  for  equation  (D-3),  the  value  of  c  must  be  calculated.   The 

appropriate  formula  is 

c2  =  nVjBCt^)  (D-7) 

The  value  of  n  is  7.15. 

Iteration  Procedure 

The  preceding  equations  permit  calculation  of  p  =  p(e,v); 

however,  an  explicit  formulation  is  not  possible,  or,  at  least,  is  very 

difficult.   Consequently,  an  iterative  approach  is  used.   Define  e  as  a 

c 

calculated  value  from  equation  (D-3).   Define  e.  as  an  input  value  at  which 

the  pressure  is  desired.  As  implied  by  Figure  D-l,  an  iterative  procedure  is 

needed  to  find  the  correct  value  for  t,  at  which  e  =  e.. 

c     1 

The  function  e   is  a  monotonically  increasing  function  in  the  temperature 
region  of  interest.   An  iteration  scheme  using  straight  lines  is  used  to  find 
t  .   For  computer  programming,  capital  letters  are  used;  hence,  in  terms  of 
capital  letters,  the  equation  for  the  straight  line  is 

F  -  F.I 
T3,J  =  E2  -  E1(T2  "  T1)  +  T1  (D"8) 

where  the  symbol  definitions  are  as  follows: 

T3,J      calculated  value  for  temperature  using  a  straight 

line  as  shown  in  Figure  D-2;  J  is  an  index  which 

indicates  the  number  of  the  iteration. 

E        desired  value  of  E  BAR;  also  identified  as  input  value  e.. 

l 

El       value  of  E  BAR  at  temperature  Tl. 

E2       value  of  E  BAR  at  temperature  T2. 
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Internal 
Energy, 


Calculated  internal  energy 
as  a  function  of  t,  e  (t) 


e. 

1 


Input  value  of 
internal  energy,  ~e, 


Temperature,  t_,   C 


Figure  D-l.   Iteration  to  determine  correct  value  of  t^ 
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STRAIGHT  LINE  -^  y 


/ 


4 

A* 

X 

E 

i 

1 

1 

1 

1 

/!     i 

=  E 


E  AS  FUNCTION 
OF  T 


Tl     T2    T3,l 
(a)  Initial  straight  line 


Tl     T2 
(b)  Calculation  of  E3,l 


SECOND  ITERATION  STRAIGHT  LINE 


.  E  AS  FUNCTION  OF  T 


Tl   T3,2   T2 


(c)  Revised  straight  line 
and  solution  for  T3,2 


(d)  Calculation  of  E3,2 


Figure  D-2.   Illustration  of  iteration  procedure, 
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Refer  to  Figure  D-2(a).   The  first  values  for  Tl  and  T2  are  inputs  and  are 
guesses.   The  program  calculates  El  and  E2  and  solves  for  T3,l  as  shown  in 
Figure  D-2(a).   The  value  for  T3,l  is  used  to  calculate  E3 .   E3  is  compared 
to  E.   If  the  difference,  E  -  E3,  is  too  large,  then  T2  =  T3,l,  Tl  =  T2, 
E2  =  E3,  and  El  =  E2.   A  new  straight  line  is  calculated  passing  through 
points  b  and  d  in  Figure  D-2(b)  and  Figure  D-2(c).   A  solution  is  obtained 
for  T3,2,  which  is  the  second  value  for  T3.   Using  T3,2,  E3  is  calculated  and 
is  once  again  compared  with  E.   If  the  quantity  ABS(E  -  E3)  is  too  large,  the 
values  for  E  and  T  are  reassigned  as  follows:   T2  =  T3,2,  Tl  =  T2,  E2  =  E3, 
and  El  =  E2. 

Once  again  a  solution  is  obtained  for  T3,3  using  a  line  passing  through 
points  d  and  e  in  Figure  D-2(d).   The  procedure  repeats  until  the  desired  accuracy 
for  ABS(E  -  E3)  is  obtained. 
Program  Description 

Figure  D-3  is  a  flow  chart  for  the  program  p  =  p(e,v),  which  has  been 
labelled  P  EV.   Additional  storage  registers  are  assigned  as  shown  in 
Table  D-l.   The  assignment  of  registers  1  to  39  is  given  in  Table  1  and 
Table  C-l.   Table  D-2  is  a  program  listing.   The  various  subroutines,  such  as 
VW,  BT,  and  H10,  must  be  in  the  memory  of  HP41CV. 

As  shown  in  the  flow  chart  of  Figure  D-3,  inputs  are  made  for  B,  E  BAR, 
J,  Tl,  and  T2.   An  input  J  =  0  must  always  be  used;  otherwise  the  program  logic 
will  be  in  disarray.   The  quantities  T  and  V  are  the  variables  used  in  subroutine 
LBL  12  to  calculate  E.   The  quantities  utilized  or  calculated  as  a  function 
of  J  are  as  follows: 

VALUE  OF  J     VALUE  OF  T     VALUE  OF  E 

0  Tl  El 

1  T2  E2 

2  T3  E3 
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INPUT  V, 

E  BAR,  J,  Tl,  T2 


LBL  12 

E  =  E(T,V) 


J=2 


LBL  11 

CALCULATE 

T3 

T  =  T3 

GTO  12 


J=2 


LBL 

14 

Tl  = 

T2 

T2  = 

T3 

El  = 

E3 

E2  = 

E3 

GTO 

11 

LBL  13 

CALCULATE 

PRESSURE 

P 


END 


Figure  D-3.   Flow  chart  for  the  program  P  EV. 


45 


Table  D-l.   Additional  Storage  Registers  for  P  EV 


Register   Symbol  Definition  Units      Programs 

40  v     Specific  volume  at  which  pressure  is     ~ 

desired  cm  /gm        12 

41  E     Internal  energy  at  which  pressure  is 

desired  kilojoule/kg  12 

42  T     Value  of  t,  as  used  in  equation  (C-5)   °C  VW,  BT 

43  -     Not  assigned 

44  v..    Value  of  specific  volume  at  p  =  0  and    „ 

temperature  t1  cm  /gm       12,  13 

45  B(t)   Parameter  in  equation  of  state 

relating  p  to  v  kilobars      12,  13 

2 

46  c.     Speed  of  sound  squared  at  temperature   2    ? 

t-  m  /sec        12 

47  h1f.   Contribution  to  enthalpy  at  zero 

pressure  between  temperatures 
tQ  and  tx 

48  z     Ratio  v  /v 

49  El    Internal  energy  evaluated  at  Tl 

50  E2    Internal  energy  evaluated  at  T2 

51  Tl    Temperature  at  point  1  on  straight 

line 

52  T2    Temperature  at  point  2  on  straight 

line 

53  J     Index  to  count  number  of  executions 

of  LBL  12  -         09 

54  E     Register  to  store  output  of  LBL  12 

which  is  E  =  E(T,v)  Joule/kg      12 


kilo. 

joule 

7kg 

12 

12 

Jouli 

e/kg 

12 

JouL 

e/kg 

12 

°C 

12 

°c 

12 

55      T3    Temperature  of  point  determined  by 
straight  line 


°C  11 


56      E3    Internal  energy  at  temperature  T3      Joule/kg      12,  14 
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Table  D-2.   Program  Listing 

01  LBLTP  EV 

26  XEQTBT 

51  1/X 

76  RCL  53 

101 

RCL  51 

126 

RCL  50 

02  TV? 

27  STO  45 

52  7.15 

77  1 

102 

- 

127 

STO  49 

03  PROMPT 

28  7.15 

53  / 

78  + 

103 

X 

128 

RCL  56 

04  STO  40 

29  X 

54  + 

79  STO  53 

104 

RCL  51 

129 

STO  50 

05  TE  BAR? 

30  RCL  44 

55  6.15 

80  GTO  11 

105 

+ 

130 

GTO  11 

06  PROMPT 

31  X 

56  1/X 

81  LBL  09 

106 

T 
T3  = 

131 

LBL  13 

07  1000 

32  1000 

57  - 

82  RCL  54 

107 

ARCL  X 

132 

7.15 

08  X 

33  / 

58  RCL  46 

83  STO  49 

108 

AVIEW 

133 

RCL  48 

09  STO  41 

34  STO  46 

59  X 

84  RCL  52 

109 

STOP 

134 

XOY 

T 
10  J? 

35  RCL  42 

60  RCL  47 

85  STO  42 

110 

STO  55 

135 

YtX 

11  PROMPT 

36  XEQTH  10 

61  + 

86  RCL  53 

111 

STO  42 

136 

1 

12  STO  53 

37  1000 

T 
62  E  = 

87  1 

112 

GTO  12 

137 

- 

13  TT1? 

38  X 

63  ARCL  X 

88  + 

113 

LBL  14 

138 

RCL  45 

14  PROMPT 

39  STO  47 

64  AVIEW 

89  STO  53 

114 

RCL  54 

139 

X 

15  STO  51 

40  RCL  44 

65  STOP 

90  CTO  12 

115 

STO  56 

140 

1.013E8 

16  TT2? 

41  RCL  40 

66  STO  54 

91  LBL  11 

116 

RCL  41 

141 

/ 

17  PROMPT 

42  / 

67  RCL  53 

92  RCL  50 

117 

- 

142 

TPRESS  = 

18  STO  52 

43  STO  48 

68  X=0? 

93  RCL  49 

118 

ABS 

143 

ARCL  X 

19  RCL  51 

44  6.15 

69  GTO  09 

94  - 

119 

50 

144 

AVIEW 

20  STO  42 

45  Y+X 

70  1 

95  RCL  41 

120 

X>Y? 

145 

END 

21  LBL  12 

46  7.15 

71  RCL  53 

96  RCL  49 

121 

GTO  13 

22  RCL  42 

47  / 

72  X>Y? 

97  - 

122 

RCL  52 

23  XEQTVW 

48  6.15 

73  GTO  14 

98  XOY 

123 

STO  51 

24  STO  44 

49  / 

74  RCL  54 

99  / 

124 

RCL  55 

25  RCL  42 

50  RCL  48 

75  STO  50 

100  RCL  52 

125 

STO  52 
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Once  J  has  acquired  a  value  of  2,  the  program  loops  repeatedly  as  follows: 
LBL  12,  TEST  ABS(E-E3),  LBL  14,  LBL  11,  LBL  12,  etc.   When  the  criterion 
for  the  value  of  ABS(E-E3)  is  attained,  the  program  departs  from  the  loop 
and  moves  to  LBL  13. 
Sample  Run 

The  program  P  EV  has  been  assigned  a  key  in  USER  mode  of  HP41CV.   Pushing 
that  key  causes  the  following  queries  to  appear: 

V?        0.75  cm3/gm 
E  BAR?     900  kilo joule/kg 
J?        0 
Tl?        50°C 
T2?        150°C 
Once  the  input  has  been  made,  the  crow's  foot  moves  across  the  display  indicating 
the  program  is  being  executed.   Next  to  appear,  after  R/S  is  pushed,  is 

E  =  377,507  (Joule/kg) 
which  is  the  value  for  El  calculated  at  Tl.  Next 

E  =  880,354  (Joule/kg) 
appears  which  is  E2.   The  situation  for  the  straight  line  is  exactly  as  shown 
in  Figure  D-2(a).   Pushing  R/S  yields 

T3  =  153.91  (=  T3,1°C) 
Pushing  R/S  causes  the  program  to  display 

E  =  902,126 
which  is  E3.   The  situation  is  that  illustrated  in  Figure  D-2(b).   The  sequence 
of  values  is  listed  below 

T3  =  153.53  (=  T3,2°C) 
E  =  899,985 
PRESS  =  30.45  (kilobars) 
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Reference  to  line  119  of  Program  Listing  indicates  that 

|E-E3|<50 

has  been  achieved. 

To  indicate  the  influence  of  the  initial  values  for  Tl  and  T2,  the 

program  was  run  again  with 

V  =  0.75  cm3/gm 

E  BAR  =  900  kilo joule/kg 

J  =  0 

Tl  =  0 

T2  =  10 

The  sequence  of  values  is  listed  below: 

E  =  135,345  (=  El) 

E  =  184,447  (=  E2) 

T3  =  155.73  (=  T3,l) 

E  =  912,390  (=  E3) 

T3  =  153.25  (=  T3,2) 

E  =  898,425 

T3  =  153.53  (=  T3,3) 

E  =  899,993 

PRESS  =  30.45 

Using  values  of  Tl  and  T2  quite  far  removed  from  final  value  of  T3,J,  the 

program  converged  to  a  solution  in  three  iterations. 

Reference 

D-l.   P.  Kutler  and  S.  R.  Chakravarthy,  "Supersonic  Flow  over  Ablated  Nosetips 
Using  an  Unsteady  Implicit  Numerical  Procedure,"  AIAA  78-213,  1978. 
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